%% Solve for the steady state of my model
% with search frictions, menu costs
% and idiosyncratic productivity
% fixed menu cost
% approximating the value of adjustment, value of not adjustment, and the
% expected value function 

%addpath ~/compecon/CEtools

soblen = 100;

% All
options.Nbell       = 100;        % Number of Bellman iterations before Newton.
options.Nnewt       = 15;       % Maximum number of Newton steps
options.tolc        = 1e-12;     % Tol on value functions
options.T_irf       = 50;       % Number of periods for IRFs
options.T           = 300;     % Number of periods for MIT shock
options.Terg        = 500;      % Number of periods to burn
options.hpisteps    = 100;

% Stationary
options.L = [];
options.itermaxL    = 5000;     % Max number of iterations to find L
options.tolL        = 1e-12;    % Tol for L
options.tolK        = 1e-5;     % Tol for equilibrium K
options.itermaxK    = 100;      % Max iterations for bisection
options.cresult = [];

options.tolD        = 1e-05;
options.itermaxp    = 25;

options.damp = 0;
options.damp_mit = .99;
options.damp_DC  = .9;

optset('bisect', 'tol', 1e-12)

%% Set globals and parameters
% % try vavra numbers
glob.n          = [30,30];   % Number of nodes for b,a
glob.nf         = [100,30];

glob.spliorder  = [1,1];    % Order of splines (Envelope Condition Method seems only robust when quadratic in k)
glob.Na         = 30;           % Number of nodes for quadrature
glob.curv       = .1;            % Amount of curvature for p grid

glob.Nsignal    = 1e6;

glob.minb = 0;
glob.maxb = 10;


%% Model parameters

% preferences
param.beta      = 0.96;               

%  rho_s    phi_L  sigma epsilon   ce_mu     d_E
% 0.60338 0.076289 2.9199   2.835 -2.3445 0.59226

% idiosyncratic shock
param.rho_s           = 0.55; % calibrate this
param.sigma_s         = 0.22;  % fix this (Clementi chooses is 0.22) 

% labor adjustment cost
param.phi_L           = 0.27084;

% Demand parameters
param.sigma           = 6.0986;              % average elasticity
param.epsilon         = 1*param.sigma;     % epsilon/sigma is superelasticity, 0 is CES 

param.omega           = 1;                  
         
param.A = 1;

% entry and exit
param.gamma = .015;     % exogenous exit rate

% fixed cost of production
param.mu_f     = -5.401;
param.sigma_f  = .9;

% sunk cost of entry 
param.ce       = exp(param.mu_f + param.sigma_f^2/2);
param.ce_noise = 1e-3;

param.M        = 10;

% dist of entrant productivity signals
param.xi     = 2;

param.nu      = 1/2; % inverse Frisch elasticity
                     % same as in clementi
                     % can calibrate aggregate shock size to hit moments

% free labor adjustment 
param.delta    = 0.1;

glob.mink       = 0;
glob.maxk       = 1e8;

%% calibrate
options.itermax_DC = 100;

P = sobolset(7);
sobolstart = (soblen-1)*sobnum + 1;
sobolend = sobolstart + soblen;

P = P(sobolstart:sobolend,:);

%    phi_L, sigma, sigma/epsilon, mu_f, sigma_f,  xi      rho_s
lb = [.001   5      .5              -6     0       1.5    .5];
ub = [.5     12     2               -3     2       4     .85];

loc = lb;
scale = ub - lb;

C = [];
M = [];

for p = 1:soblen
    calib = P(p,:);
    calib = loc + scale.*calib
    [mom] = compute_moments_for_calibration(calib, param, glob, options);

    if flag
       continue 
    end
    
    dlmwrite(['results/calib', num2str(sobnum), '.csv'],calib,'delimiter', ',','-append');
    
    f = fields(mom);
    mom_out = nan(length(f),1);
    for i = 1:length(f)
        mom_out(i) = mom.(f{i});
    end
    
    dlmwrite(['results/moments', num2str(sobnum), '.csv'],mom_out','delimiter',',','-append');

end

